REM ------ Q2iteration ---------Iteration of a rational function ---------
REM
REM           Q2(z)=z/(z^2-az+a),   z=x+iy, a=ar+i*ai
REM  
REM         Z0（白の点）を与えての反復計算，Zi (i=1,2,3,...) は複素平面上に
REM         プロットする。カラーは3色使う：緑，黄，赤
REM
REM--------------------------- Joe  5.25 '19------------------------------   
DIM co(3),it(0 TO 3)
READ co(1),co(2),co(3)
DATA 3,6,4
READ txl,txu,tyl,tyu
DATA .1,.9,.1,.9
LET it(0)=0
LET EPS=.0001
LET IS=0
LET ia=0
INPUT PROMPT "パラメタaのinput (ar,ai)=":ar,ai
print"複素平面上の四角形のinput"
INPUT PROMPT "xの範囲 xl,xu=":xl,xu
INPUT PROMPT "yの範囲 yl,yu=":yl,yu
REM ---- 画像パラメタ設定  ----
LET wx=xu-xl
LET wy=yu-yl
LET twx=txu-txl
LET twy=tyu-tyl
LET a1=twx/wx
LET a2=twy/wy
LET b1=(xu*txl-xl*txu)/wx
LET b2=(yu*tyl-yl*tyu)/wy
SET AREA COLOR 1
PLOT AREA:txl,tyl;txu,tyl;txu,tyu;txl,tyu;txl,tyl
SET LINE COLOR 7
IF xl*xu<0 THEN PLOT LINES:b1,tyl;b1,tyu
IF yu*yl<0 THEN PLOT LINES:txl,b2;txu,b2
REM --------
INPUT PROMPT "初期値 x0,y0=":x0,y0
INPUT PROMPT "iterationの回数 tm=":tm
PLOT TEXT, AT .05,.96:"Q2iteration, a="
PLOT TEXT, AT .25,.96,USING"###.###### ###.######":ar,ai
PLOT TEXT, AT .05,.92:"[xl,xu]*[yl,yu]="
PLOT TEXT, AT .3,.92,USING"##.##### ##.#####  ##.##### ##.#####":xl,xu,yl,yu
PLOT TEXT, AT .7,.04:"反復回数 tm="
PLOT TEXT, AT .84,.04,USING"#####":tm
LET it(3)=tm
LET it(1)=INT(tm/3+.5)
LET it(2)=INT(2*tm/3+.5)
SET POINT STYLE 3
SET POINT COLOR 0
LET tx0=a1*x0+b1
LET ty0=a2*y0+b2
PLOT  POINTS: tx0,ty0
REM
REM     計算開始
REM
LET x=x0
LET y=y0
FOR i=1 TO tm
   LET d1=x*x-y*y-ar*x+ai*y+ar
   LET d2=2*x*y-ai*x-ar*y+ai
   LET d=d1*d1+d2*d2
   IF d<EPS THEN LET IS=IS+1
   LET x1=(x*d1+y*d2)/d
   LET y1=(y*d1-x*d2)/d
   LET d3=x1*x1+y1*y1
   IF i>100 AND d3>16 THEN LET ia=ia+1
   PRINT i;")";x1;", ";y1
   IF it(0)<i AND  i<=it(1) THEN  SET POINT COLOR co(1)
   IF it(1)<i AND  i<=it(2) THEN  SET POINT COLOR co(2)
   IF it(2)<i AND  i<=it(3) THEN  SET POINT COLOR co(3)
   REM
   LET tx=a1*x1+b1
   LET ty=a2*y1+b2
   PLOT  POINTS: tx,ty
   LET x=x1
   LET y=y1
NEXT i
PRINT"分母ほとんどゼロの回数=";IS
PRINT"n>100 で絶対値4以上の回数=";ia 
REM
REM ---------------------------
INPUT PROMPT"ファイル名，コメント書くか? y/n ":q$
IF q$="n" OR q$="n" THEN GOTO 99
PRINT "ファイル名などinputせよ(全角20文字以下) (例）Q2it-01,Q2it-02,.."
LINE INPUT PROMPT "ファイル名など:":na$
PLOT TEXT, AT .04,.04:na$
99 END